Vector solitons in nearly-one-dimensional Bose-Einstein condensates 
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We derive a system of nonpolynomial Schrodinger equations (NPSEs) for one-dimensional wave 
functions of two components in a binary self-attractive Bose-Einstein condensate loaded in a cigar- 
shaped trap. The system is obtained by means of the variational approximation, starting from 
the coupled 3D Gross-Pitaevskii equations and assuming, as usual, the factorization of 3D wave 
functions. The system can be obtained in a tractable form under a natural condition of symmetry 
between the two species. A family of vector (two-component) soliton solutions is constructed. 
Collisions between orthogonal solitons (ones belonging to the different components) are investigated 
by means of simulations. The collisions are essentially inelastic. They result in strong excitation 
of intrinsic vibrations in the solitons, and create a small orthogonal component ( "shadow" ) in each 
colliding soliton. The collision may initiate collapse, which depends on the mass and velocities of 
the solitons. 

PACS numbers: 03.75.Ss,03.75.Hh,64.75.+g 



T3 

C 

O 
o 



O 

-i— » 



i 

C 

o 
o 



X 



I. INTRODUCTION 

Bose-Einstein condensates (BECs) with attractive in- 
teractions between atoms can form stable wave packets in 
nearly one-dimensional (ID) "cigar-shaped" traps, which 
provide for tight confinement in two transverse direc- 
tions, while leaving the condensate almost free along the 
longitudinal axis. This trapping geometry made it possi- 
ble to create stable bright solitons and trains of such 
solitons in the 7 Li condensate, in which the interaction 
between atoms was made weakly attractive by means of 
the Feshbach-resonance technique. In the 85 Rb conden- 
sate trapped under similar conditions, stronger attraction 
between atoms leads to controllable collapse and creation 
of nearly 3D solitons 0] • 

This experimentally relevant situation is described by 
effective ID equations which may be derived from the 
full 3D Gross-Pitaevskii equation (GPE) under various 
conditions and by means of different approximations Q- 
@. In some cases, the deviation of the effective equation 
from a straightforward ID variant of the GPE amounts to 
keeping an extra self-attractive quintic term in the equa- 
tion, which may be sufficient to essentially alter proper- 
ties of the corresponding solitons jy, fej, ■ A more con- 
sistent derivation, that starts with the factorization of the 
3D wave function into the product of a transverse one (it 
actually represents the ground state of the 2D harmonic 
oscillator) and arbitrary slowly varying longitudinal (one- 
dimensional) wave function, and then uses the variational 
approximation ^l]], leads to a more sophisticated but 
also more accurate nonpolynomial Schrodinger equation 
(NPSE) for the longitudinal wave function 0, H2- The 
above-mentioned simplified equation featuring the com- 
bination of cubic and quintic terms can be then obtained 
by an expansion of the NPSE for the case of a relatively 
weak nonlinearity pJJ. The ratio of the coefficients in 
front of the cubic and quintic terms in the model derived 



in Ref. |6j is not the same as follows from the expansion 
of the NPSE, which is explained by a coarser character 
of the approximation used in that work (the approxima- 
tion did not allow a deviation of the nonlinearity different 
from the cubic-quintic form). 

A physically significant generalization of the above- 
mentioned equations is a system of two nonlinearly cou- 
pled equations for a binary BEC, which can be created in 
the ex peri ment by means of the sympathetic- cooling tech- 
nique [l3|. Accordingly, a relevant problem is to derive 
a system of effective ID equations for a mixture of two 
BEC species in the cigar-shaped trap, starting from the 
two coupled GPEs in the 3D space. In this work, we aim 
to derive such a system in the form of coupled NPSEs, 
using a gene ralized version of the method elaborated in 
Refs. Ull El- 

The paper is organized as follows. The derivation of 
the coupled NPSE system, which is based on the varia- 
tional approximation, is presented in Section II. In the 
most general case, it leads to a cumbersome system. 
However, we demonstrate that, under a natural condi- 
tion of a symmetry between the two species, the equa- 
tions may be reduced a tractable closed system of two 
NPSEs for longitudinal wave functions of the two com- 
ponents. Then, in Section III, we consider solutions for 
vector solitons (i.e., two-component ones) generated by 
this system; the solutions are found in an implicit ana- 
lytical form up to a point where they cease to exist due 
to collapse. 

A natural application of the thus derived NPSE system 
is to consider collisions between two orthogonal solitons, 
which belong to the two different components. This anal- 
ysis, based on numerical simulations, is reported in Sec- 
tion IV. The collisions arc inelastic, which is manifested 
in the excitation of intrinsic oscillations in the solitons 
after the collision, and generation of small "shadows" in 
them (each soliton captures and keeps a small share of 
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atoms from the other species). The strongest manifes- 
tation of the inelasticity, as we demonstrate in Section 
IV, is a possibility to initiate collapse by the collision be- 
tween two orthogonal solitons (which depends on their 
relative velocity). The paper is concluded by Section V. 



II. COUPLED NONPOLYNOMIAL 
SCHRODINGER EQUATIONS 

The system of 3D GPEs for a dilute binary condensate, 
confined in the transverse direction by a strong harmonic 
potential with frequency u± and in the axial direction by 
a generic weak potential V(z), can be derived from the 
Lagrangian density, 



L = 



fc=l,2 



2 V 



y 2 ) 



-V(z) - irg k \fa\ 2 ] fa ~ 2n g 12 \fa\ 2 \fa\ 



(1) 



Here 'i/'fc (r, i) is the macroscopic wave function of the k- 
th species, which is subjected to the normalization con- 
dition, 



where real <r k (z,t) and complex f k (z,t) are dy- 
namical fields, the latter ones obeying normalization 
J-°° \fk(z)\ 2 dz = N k , as it follows from Eqs. J2J and 

Inserting this ansatz in Lagrangian density |JT)|. per- 
forming the integration in the transverse plane, and ne- 
glecting derivatives of cr(z, t) (for the same reasons as in 
Refs. 13 )> one can derive the following effective La- 
grangian: 



k=l,2 



1 



id t +-di--\ 



1 



4 



-V(z) 



Iff* 
2a 



/fc-2 



■9l2|/l| 2 |/ 2 P 



This Lagrangian gives rise to a system of four Eulcr- 
Lagrange equations, obtained by varying L with respect 
to /£ and a k : 



idtfk 



1 



d\ + v(z) + 
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\fa(x,y,z)\ dxdydz = N k , 



(2) 



where Nk is the number of atoms in the fc-th species, and 



9k = 2a k /a^, g 12 = 2a 12 ja\_ 



(3) 



are strengths of the intra- and inter-species interactions, 
where a k and a± 2 are the scattering lengths, and a± — 
yh] (muj±) the transverse harmonic-confinement length. 
Here we consider the binary condensate with attraction 
between atoms, which implies that both 0,1,2 and 012 are 
negative. In the Lagrangian density, lengths, time, and 
energy are written in units a±, wj 1 , and hco±, respec- 
tively. 

The ordinary variational procedure applied to Eq. l|Tjl 
gives rise to the coupled 3D GPEs, 
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(7) 



with k = 1,2 in both equations 10 and J7J). This is 
a full system of the coupled NPSEs describing the two- 
component nearly-lD BEC. 

Note that the ansatz (JSJ is relevant when the transverse 
confinement size is much smaller than a characteristic ax- 
ial length of a structure (in particular, soliton) to be ob- 
tained as a solution of the axial equation. Physically, this 
means that the quantum pressure in the transverse direc- 
tion is much stronger than the nonlinear self-attraction 
in the condensate. 



+2Tr 5A ,# fe | 2 + 27r 3l 2^ 3 -fe| 2 ] fa, k = 1,2. (4) 

Our objective here is to derive a system of effective ID 
NPSEs, following the lines of the derivation of the NPSE 
for the single-component condensate developed in Ref. 
[its generalization for an axially nonuniform trapping 
potential, with u± = uj±(z), was reported in Ref. 0]- 
Using the cylindrical coordinates (r, 8) in the transverse 
plane (x, y), we adopt the usual ansatz for the wave func- 
tions strongly localized in this plane, and weakly confined 
in the axial direction, z: 



fa(r,z,t) 



1 



/7TCTfc(z,t) 



exp 



2a k (z,t) 



2 r fk(z,t) , 



(5) 



III. VECTOR BRIGHT SOLITONS 

In the subsequent analysis of the coupled NPSEs, we 
focus on the most natural symmetric case, when the (neg- 
ative) effective nonlinearity coefficients accounting for the 
intra- and inter-species interactions are equal, namely, 



512 = 9i = 92 



(8) 



As follows from Eqs. ©, this relation takes place, in 
particular, when the all scattering lengths are equal. In 
the symmetric case, we assume that numbers of atoms 
in the two species are equal too. We note that, for the 
self- repulsive binary BEC, with g > (recall here we 
are going to consider the case of g < 0), Eq. © may 



3 



pose a formal problem, as it precisely corresponds to the 
miscibility-immiscibility threshold in the infinite system. 
Nevertheless, the problem does not really take place, as 
the pressure exerted by the external potential shifts the 
equilibrium towards the miscibility (see, e.g., Ref. [l5jp. 
hence the case corresponding to relation (JHJ) in the re- 
pulsive binary BEC is not going to be a degenerate (i.e., 
structurally unstable) one. 

If constraint JSJ) holds, Eqs. (J7| take the form 



vi = l + g\f k \ 2 + Hh-k\ 



and admit an exact symmetric solution: 



(9) 



(10) 



Of course, there remains a question if some additional 
asymmetric solutions to Eqs. © may also exist. One 
may assume that an asymmetric solution, if any, branches 
off from the symmetric one through a bifurcation. Then, 
close to the bifurcation point, one will have a\ 2 — 
<7q (1 + #1,2)1 with some infinitesimal 5\ 7^ 82 [o~o is the 
symmetric solution given by Eq. HlOjl ] . Substituting this 
in Eqs. JJjJ and linearizing them in <5i and 82, one arrives 
at a system 



24 
F k 



F k (S k 
9\fk\ 2 /4 



3-fcJ 



1,2 



(11) 
(12) 



The resolvability condition for linear system (|ll(l (equat- 
ing its determinant to zero) takes the following form, af- 
ter simple calculations: F\ + F2 = 2. However, this con- 
dition cannot hold for the attractive BEC, with g < 0, 
because expressions F\ and F 2 , as given by Eq. fl^|) . are 
negative in this case. This means the bifurcation giving 
rise to asymmetric solutions is impossible in the attrac- 
tive binary condensate [provided that constraint JHJ) is 
valid] , which substantiates the use of symmetric solution 

m- 

The substitution of Eq. l(TU|) in Eqs. JBJ leads to 
closed-form equations for the complex amplitude func- 
tions, /1 and f 2 , 



.dfk 
''~dt 



1 d 2 
-2^ + V ^ 



l/il 2 



l/ 2 | 2 



yi+.9(i/ii 2 + i/ 2 i 2 ) 



+ V / l + 5(l/i| 2 + l/2| 2 ) 



(13) 

Equations l|13|) reduce to the familiar integrable Man- 
akov's system (MS) [T^ . 



idtfk = 



--dl + V{z)+g(\h\ 2 + \f2\ 2 ) 



.fk 



(14) 



if J((|/i 1 2 + I/2 1 2 ) *C 1. Only under this condition 
the system may be considered as truly one-dimensional, 



and only in this limit it is integrable. Nevertheless, 
in the general case Eqs. (fH?)) share the "isotopic in- 
variance" with the Manakov's system: the nonlinear- 
ity appears solely through the invariant combination, 
|/i| 2 + |/ 2 | 2 - Due to this fact, Eqs. (|13|) conserve an 
additional dynamical invariant ("isotopic spin"), S = 

jl^ [A (s)/! (z) + fi(z)fo(z)] dz -, with asterisk standing 
for the complex conjugate. 

In the case of attraction, g < 0, vector (two- 
component) bright solitons arc looked for as fk = 
exp (— ifikt) $k(z), where <E>i and $2 are real localized 
functions obeying the following coupled equations: 



V(z)$ k 
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l+.g($ 2 + a>2) 



4> 



k ■ 



(15) 

Due to the isotopic invariance, the solitons with equal 
chemical potentials of their components are tantamount 
to the single-component (scalar) one, with $2=0 and 
$1 = $(z) being a solution of a single equation, 
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(16) 



If a soliton solution to Eq. (|16fl is found, the correspond- 
ing vector soliton can be constructed in an obvious way, 



A(M) 
A(M) 
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exp (—i/it) 
exp (—ijit) 



(17) 



with an arbitrary "isotopic angle", < 9 < tt/2. More 
general vector solitons, with different chemical potentials 
in their components, are possible two. However, in the 
symmetric case that we are dealing with here, it is ob- 
vious that the vector solitons with unequal chemical po- 
tentials cannot realize the ground state, therefore they 
are considered here. 

For V(z) = and g < 0, a family of soliton solutions 
to Eq. (fTfj|l was constructed in Refs. I n this case, 

setting $(z) = \/Ntf>(z), with Ni = N 2 = N, and 



1=N\9\ , 



(18) 



the field <j)(z) and the chemical potential /i are given by 
implicit formulas, 



1 
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l + /x 
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(20) 



and the wave function satisfies the normalization con- 
dition fj 1-00 (</>(z)) 2 dz = 1. The family is then charac- 
terized by the dependence of 7 on /1. The inverse of Eq. 
(|20|l demonstrates that, in terms of the /z(7) dependence, 
there are two branches of the soliton family, but only one 
of them, that satisfies condition dfi/dj < (which is 
nothing else but the known Vakhitov-Kolokolov stability 
criterion [2(J in the present notation), is stable. In ad- 
dition, there is a critical nonlinearity strength, 7 C = 4/3 
(which corresponds to fj, = 1/2), above which the solu- 
tion does not exist, because of the longitudinal collapse 
0, ^| , which is a manifestation of the 3D collapse pos- 
sible in the underlying system of GPEs, Eqs. In the 
limit of weak nonlinearity 7^0, Eq. (|19|) reduces to the 
ordinary soliton waveform, </>(z) = sech (72/2) 
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FIG. 1: Peak density np of the two colliding solitons as a func- 
tion of time (t), in the nonintegrable system of NPSEs, Eqs. 
(1131 . and in the integrable Manakov's system (alias NLSEs), 
Eqs. 1141 . In both cases, the initial velocity is v = 0.6. 



IV. COLLISIONS BETWEEN SOLITONS 

A straightforward application of the system of NPSEs 
()13JI is to study collisions between two orthogonal soli- 
tons, taken in the form of Eq. I|17|l . with equal values 
of [i and isotopic angles 6 — and 9 — tt/2. Using the 
Galilean invariance of the equations, the velocities of the 
two solitons are taken to be ±i>. The corresponding ini- 
tial condition, at t = 0, is 

j <p[ 0) (z)\ _ / \ , 211 

with large initial separation zq. To determine the time- 
evolution of the fields <pk{z,t), k = 1,2, we solved 
both NLSEs and NPSEs numerically, by using a two- 
component extension of a well-tested finite-difference 
code based on the Crank-Nicolson predictor-corrector al- 
gorithm [2l|] . In the MS [alias the NLSE system, Eqs. 
(|14fl ], which is integrable, collisions are always elastic. 
However, since NPSEs <|13fl are not integrable, collisions 
described by these equations are expected to be inelas- 
tic. This expectation is borne out by Fig. 1, where we 
compare the collision outcomes in the MS and NPSEs for 
identical sets of parameters. The figure shows the peak 
densities, np, of both colliding solitons (which are equal, 
due to the symmetry of the configuration) , as a func- 
tion of time. The outcome does not depend on the initial 
separation zq between the solitons in Eq. (|2ip. provided 
that it is large enough (we took zq = 200) . After the col- 
lision the MS solitons remain undisturbed (dashed lines), 
while their NPSE counterparts come out from the colli- 
sion with excited intrinsic oscillations (solid lines). This 
result not only shows that the collision in the NPSEs is 
inelastic, but also suggests that the solitons supported by 
this system, i.e., ones given by Eqs. ((T7|l . ijT§|) . and (|2U|l . 
unlike their counterparts in the integrable MS, feature 
an intrinsic mode, with a well-defined eigenfrequency, u>. 
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FIG. 2: Peak density np of the two colliding solitons as 
a function of time, found from numerical integration of the 
NPSEs, Eqs. I13L for different values of interaction strength 
7. The initial velocity is v — 0.6. 



In fact, this mode was predicted in Ref. |12(, by means 
of the variational approximation. It was shown that u> 
vanishes at 7 — > 0, and it attains a maximum close to 
the above-mentioned collapse threshold, 7 = 7 C = 4/3. 

As shown in Fig. 2, the amplitude and frequency (uj) 
of the oscillations excited by the collisions of solitons in 
the NPSE system grow with interaction strength 7. On 
the contrary to that, the simulations demonstrate that 
the amplitude and frequency of the intrinsic oscillations 
do not depend on initial velocity v (see Fig. 5 below). 
The independence of u> on v is quite natural, as the fre- 
quency is determined solely by the internal structure of 
the soliton. 

In Fig. 3 we compare the intrinsic frequency, to, as 
found from the direct simulations of the NPSEs, Eqs. 
(|13|l . to the frequency calculated by means of the varia- 
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FIG. 3: Frequency ui of the intrinsic oscillations of the two 
solitons as a function of interaction strength 7. Stars: ui found 
from the numerical solution of NPSEs, Eqs. 11311 . as the fre- 
quency of intrinsic oscillations excited by the collision between 
two solitons. Solid line: a result of the variational approxi- 
mation from Ref. [T^|. The collapse point, found from simu- 
lations of Eqs. Ill .'it . is indicated by the dashed vertical line. 
The initial velocity is v = 0.8. 



FIG. 4: The trapping effect in the collision of two solitons: 
small parts of fields (j>2 {z, t) and (j>i (z, t) remain bound, respec- 
tively, in the first and second soliton after the collision. Top, 
central, and bottom panels display the density configurations, 
n ia {z) ee |0i,2 (^) | 2 , at t = 100, 180, and 240. The interaction 
strength is 7 = 1.2, and the initial velocity is v = 0.6. 
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tional approach of Ref. 12] . The figure shows that the 
variational approximation somewhat overestimates both 
the critical strength of longitudinal collapse, 7 C , and fre- 
quency uj. For both quantities, the relative error is about 
15%. 

Collisions between solitons in NPSEs give rise to an- 
other noteworthy effect: after the collision, a small part 
of the field, 4>i , which originally belonged to the first soli- 
ton remains trapped in the second soliton, and vice versa, 
see Fig. 4. The effect is visible only for 7 > 1, in the 
velocity interval of 0.4 < v < 0.7 (outside this parameter 
range, the effects takes place too but is very weak). Note 
that a similar effect ("shadow formation") was observed 
in the model describing the interaction of two polariza- 
tions of light in a nonlinear optical fiber, which was based 
on a nonintegrable system of two NLSEs with the cubic 
nonlinearity, see Refs. [T|| H9( and references therein. 
An explanation of the trapping effect is based on the fact 
that the soliton in each field <fik(z,t) (fc = 1,2) supports 
not only the above-mentioned mode of intrinsic vibra- 
tions in the same field, but also an external eigenmode 
of small perturbations in the mate field, (fe-fc(z, t). The 
latter mode is excited as a result of the collision |18|. 

The collisions feature a trend to be more inelastic at 
smaller velocities, as illustrated by Fig. 5, which shows 
the maximum and minimum values of the peak density, 
rip ^ and np^ , in the oscillating solitons emerging from 
the collision. The figure shows that, while n-p 1 ^ does not 

depend on initial velocity v, is smaller at smaller 

velocities, which implies stronger inelasticity. The results 
are shown in Fig. 5 for 7 = 1, and similar trends are 
found for 7 < 1. For completeness, in Fig. 5 we also plot 
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FIG. 5: Dependence on the collision velocity v of the maxi- 
mum (triangles) and minimum (rhombuses) values of the peak 
density, rip and n^™\ in solitons disturbed by the collision. 
The respective frequency of the intrinsic mode, uj, is shown 
by stars. The interaction strength is 7 = 1. 



frequency uj of the intrinsic-mode excited by the collision, 
which confirms that uj does not depend on v. 

An interesting issue is whether the collision may re- 
sult in collapse. As follows from Eqs. I|13fl . the collapse 
happens when condition |<?|(|/i| 2 + I/2I 2 ) = 1 takes place 
at some point. If the first maximum of the peak den- 
sity is achieved when the centers of the colliding soli- 
tons nearly coincide, this condition can be estimated as 
rip 1 ' ~ 1/(27). While Fig. 5 shows that rip ' does not 
depend on collision velocity v at 7 < 1, it depends on v 
at 7 > 1, and the collapse can thus been reached by in- 
creasing v. In particular, Fig. 6 shows that, for 7 = 1.2, 
the maximum value of peak density, np(t), grows with 
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FIG. 6: Peak density np of the two colliding solitons as a 
function of time for different values of initial velocity v. The 
collapse is induced by the collision at v = 0.9 (the corre- 
sponding curve shoots up vertically at the collapse moment, 
t = t c fH 110). The interaction strength is 7 = 1.2, and the 
initial separation is zq = 200. 



tions (NPSEs) for longitudinal wave functions of two 
components in a binary BEC, in the case of attraction 
between the atoms. The system was derived by means of 
the variational approximation, starting from the coupled 
3D Gross-Pitaevskii equations for the two species and 
assuming (as usual) the factorization of 3D wave func- 
tions into products of the strongly confined transverse 
and slowly varying longitudinal ones. The system was 
cast in a tractable form under a natural symmetry con- 
straint. Then, a family of two-component (vector) soliton 
solutions was obtained, and collisions between orthogo- 
nal solitons (each belonging to one component only) were 
studied in detail by dint of systematic numerical simula- 
tions. It was found that the collisions are inelastic. They 
lead to strong excitation of intrinsic oscillations in the 
solitons emerging from the collision, and to formation of 
a small orthogonal component ("shadow") in each soli- 
ton. Eventually, the collision may initiate collapse of the 
solitons, depending on their mass and velocities. 
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v, and the collapse takes place at v = 0.9. Note that 
collapse induced by the collision between two solitons in 
a single NPSE was reported in Ref. but in that case 
the onset of the collapse did not depend on the initial 
velocity. 

V. CONCLUSIONS 

In this work, we have derived a system of one- 
dimensional coupled nonpolynomial Schrodinger equa- 
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